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ABSTRACT 

Context. Gravitational fragmentation has been proposed as a mechanism for the formation of giant planets in close orbits around solar-type 
stars. However, it is debatable whether this mechanism can function in the inner regions (R < 40 AU) of real discs. 
Aims. We investigate the thermodynamics of the inner regions of discs and their propensity to fragment. 

Methods. We use a newly developed method for treating the energy equation and the equation of state, which accounts for radiative transfer 
effects in SPH simulations of circumstellar discs. The different chemical and internal states of hydrogen and the properties of dust at different 
densities and temperatures (ice coated dust grains at low temperatures, ice melting, dust sublimation) are all taken into account by the new 
method. 

Results. We present radiative hydrodynamic simulations of the inner regions of massive circumstellar discs and examine two cases: (i) a disc 
irradiated by a cool background radiation field (T BGR ~ 10 K) and (ii) a disc heated by radiation from its central star (T BGR ~ ^OOK^/AU] -1 ). 
In neither case does the disc fragment: in the former because it cannot cool fast enough and in the latter because it is not gravitationally 
unstable. Our results (a) corroborate previous numerical results using different treatments for the hydrodynamics and the radiative transfer, and 
(b) confirm our own earlier analytic predictions. 

Conclusions. Disc fragmentation is unlikely to be able to produce giant planets around solar-type stars at radii < 40 AU. 

Key words. Accretion, accretion discs - Hydrodynamics - Instabilities - Radiative Transfer - Planetary systems: formation - Planetary 
systems: protoplanetary discs 

1- Introduction 

A large number of exoplanets (~ 265) have been discovered 
around nearby stars in the last 12 years, mainly with the ra- 
dial velocity method (Udry & Santos 2007; Extrasolar plan- 
ets Encyclopaedia at |http : //ex oplanet . eu). Most of these 
exoplanets have properties that are different from those in the 
Solar System. In particular, giant planets are found on very 
close orbits of apparently random eccentricity. 

Two main theories have been proposed for the formation 
of giant planets, (i) core accretion, and (ii) gravitational frag- 
mentation. In the core accretion scenario giant planets form by 
coagulation of planetesimals (e.g. Safronov 1969; Goldreich & 
Ward 1973; Pollack et al. 1996). Once a solid body of around 
10 M is reached, it quickly accretes a large gaseous envelope. 
One of the main problems with this theory is that the timescale 
for planet formation is long. The theory predicts that planets 
can form within a few million years, but observations suggest 
that circumstellar discs may not be that long lived (Haisch, 
Lada & Lada 2001 ; Cieza et al. 2007). In the gravitational frag- 
mentation scenario, giant planets form by gravitational insta- 



bility in massive discs (e.g. Kuiper 1951; Cameron 1978; Boss 
1997; Durisen et. al 2007). The main advantage of the gravi- 
tational fragmentation theory is that planets condense out on a 
dynamical timescale, i.e. ^ 10 3 years. 

There are two conditions that must be fulfilled for discs 
to fragment gravitationally. (i) The disc must be gravitation- 
ally unstable, i.e. massive enough so that gravity can overcome 
thermal pressure and centrifugal support (Toomre 1964), 

where c is the isothermal sound speed, k is the epicyclic fre- 
quency, X is the surface density, and R is the distance from 
the axis of rotation, (ii) A proto-fragment must be able to cool 
fast enough for the energy delivered by compression to be ra- 
diated away. Theory and simulations (Gammie 2001; Johnson 
& Gammie 2003; Rice et al. 2003, 2005; Mayer et al. 2004; 
Mejia et al. 2005) indicate that the cooling must happen on a 
dynamical time-scale, 

W < C(r)r 0RB , o.5 £C(y)£ 2.0, (2) 
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where t 0RB is the local orbital period, and y the adiabatic expo- 
nent. 

However, it is uncertain whether real discs actually satisfy 
the above conditions. Boss (2004) and Mayer et al. (2007) sug- 
gest that convection provides the necessary cooling, whereas 
Johnson & Gammie (2003), Boley et al. (2006) and Nelson 
(2006) find that the cooling is too slow, and hence fragmenta- 
tion is not possible. The latter point of view is supported by an- 
alytic studies which indicate that convection cannot provide the 
required cooling (Whitworth et al. 2007; Rafikov 2007), and 
that discs cannot cool fast enough, at least close (R < 50 AU) 
to the central star (Rafikov 2005; Matzner & Levin 2005; 
Whitworth & Stamatellos 2006). 

One reason why hydrodynamic simulations have produced 
contradictory results concerning whether disc fragmentation 
can produce giant planets close to the central star is the dif- 
ferent treatments of radiative transfer used. We have recently 
developed an efficient scheme to capture the thermal and radia- 
tive effects when protostellar gas fragments (see Section 2 and 
Stamatellos et al. 2007). Thus, we are able to perform radiative 
SPH simulations of the inner disc regions and to include the 
effects of the equation of state consistently, i.e. by solving the 
relevant Saha equations and taking into account the resulting 
continuous and differentiable changes in the mean molecular 
weight and the internal energy. Moreover our radiative scheme 
allows us to include irradiation of the disc by a background ra- 
diation field. Our simulations suggest that it is very difficult for 
planets to form by gravitational instability in the inner regions 
of a massive circumstellar discQ 

The structure of the paper is as follows. In Section 2 we de- 
scribe our radiative hydrodynamic code. In Section 3 we define 
the initial conditions. In Section 4 we present and discuss the 
simulations. In Section 5 we summarise the results and their 
implications for the possibility of forming giant planets close 
to a star by gravitational fragmentation. 

2. Numerical method 

For the hydrodynamics we use the SPH code dragon (Goodwin 
et al. 2004), which invokes an octal tree (to compute grav- 
ity and find neighbours), adaptive smoothing lengths, multiple 
particle timesteps, and a second-order Runge-Kutta integration 
scheme. The code uses time-dependent viscosity with param- 
eters a* = 0.1, p - 2a (Morris & Monaghan 1997) and a 
Balsara switch (Balsara 1995), so as to reduce artificial shear 
viscosity (Artymowicz & Lubow 1994; Lodato & Rice 2004; 
Rice, Lodato, & Armitage 2005). 

The energy equation is treated with the method of 
Stamatellos et al. (2007). This method uses the density and the 
gravitational potential of each SPH particle to define a pseudo- 
cloud around each particle, through which the particle cools 
and heats. The mass- weighted mean column-density £/, and the 
Rosseland-mean opacity K R (pi,Ti), averaged over every pos- 
sible position of the particle inside its pseudo-cloud, are then 



used to calculate the net radiative heating rate for the particle, 
according to 



1 For simulations of gravitational fragmentation in the outer regions 
of a massive circumstellar disc, the reader is referred to Stamatellos, 
Hubber & Whitworth (2007). 



dui 
dt 



4 (t (T 4 - T 4 ) 

SB v BGR _ j * 

2^K R (p h Ti)^K p - l (p h Ti) 



(3) 



here, the positive term on the right hand side represents heat- 
ing by the background radiation field, and ensures that the gas 
and dust cannot cool radiatively below the background radia- 
tion temperature T BGR . K p (pi, T t ) is the Planck-mean opacity, 
and <r SB the Stefan-Boltzmann constant. 

The method takes into account compressional heating, vis- 
cous heating, radiative heating by the background, and ra- 
diative cooling. It performs well, in both the optically thin 
and optically thick regimes, and has been extensively tested 
(Stamatellos et al. 2007). In particular it reproduces the de- 
tailed 3D results of Masunaga & Inutsuka (2000), Boss & 
Bodenheimer (1979), Boss & Myhill (1992), Whitehouse & 
Bate (2006), and also the analytic test of Spiegel (1957). It is ef- 
ficient and easy to implement in particle- and grid-based codes. 
Because the diffusion approximation is applied here globally, 
the method does not capture in detail the exchange of heat be- 
tween neighbouring fluid elements. Our simulations also have 
insufficient resolution to capture convection. 

The gas is assumed to be a mixture of hydrogen and he- 
lium. We use an equation of state (Black & Bodenheimer 1975; 
Masunaga et al. 1998; Boley et al. 2007a) that accounts (i) for 
the rotational and vibrational degrees of freedom of molecular 
hydrogen, and (ii) for the different chemical states of hydro- 
gen and helium. However, we note that in the simulations pre- 
sented here the temperature never becomes hot enough for sig- 
nificant dissociation of H 2 , and consequently the mean molec- 
ular weight is approximately constant. We assume that ortho- 
and para-hydrogen are in equilibrium. 

For the dust and gas opacity we use the parameterization 
by Bell & Lin (1994), ic(p,T) = k p a T b , where k , a, b 
are constants that depend on the species and the physical pro- 
cesses contributing to the opacity at each p and T. The opac- 
ity changes due to ice mantle melting, the sublimation of dust, 
molecular and H" contributions, are all taken into account. 

3. Disc initial conditions 

We simulate a 0.07 Mq disc around a 0.5 Mq star. This is a rel- 
atively massive disc, but such discs have been observed, for ex- 
ample in the Orion Nebula cluster (Eisner & Carpenter 2006). 
The disc initially extends from 2 to 40 AU. Its initial surface 
density and temperature are 

/ R \~V 2 



T(R) = 



(4) 



1/2 



(5) 



where X = 1.8 x 1O" 4 M AU" 2 and T = 1200 K are the 
surface density and temperature at R 1 AU, and Too = 10 K 
is the asymptotic temperature far from the central star. These 
are typical disc profiles and have been chosen so as to match 
approximately the initial density and temperature profiles used 
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in the simulations of Boley et al. (2006) and Cai et al. (2007). 
We assume that the initial disc temperature is independent of 
distance from the disc midplane. 

To calculate the initial disc thickness, Zo(R), we balance the 
vertical gravitational acceleration due to the star and the under- 
lying disc, against the hydrostatic acceleration, 



GM+ zo(R) 



R 2 

zo(R) 



R 



+ nGL(R) 



ttI i (R)R 3 
2M+ 



c\R) 
zq(RY 

7TZ(R)R 2 

2M+ 



R 5 



-c\R) 



1/2 



(6) 



(7) 



where c(R) is the local isothermal sound speed. 

For the initial vertical structure of the disc, we adopt a si- 
nusoidal profile 



p(R,z)=p(R, 0)cos 



nz 



2z (R) 



\z\<z (R). 



(8) 



This profile resembles a Gaussian profile (e.g. Frank, King & 



Raine 1992). Setting X(R) 
7rSo 



p(R,z) 



4R l/2 zo(R) 



cos 



■zo(R) 
zo(R) 

71 Z 



2zo(R) 



p(R, z)dz , we obtain 

\z\<zo(R). (9) 



The choice of the initial density and temperature pro- 
files is not critical, since the disc quickly relaxes to a quasi- 
equilibrium state. 

4. Simulations 

We perform two simulations, one with a uniform blackbody 
background radiation field having temperature Tbgr = 10 K 
(cf. Boley et al. 2006; Cai et al. 2007), and one which takes 
account of radiation from the central star (Tbgr ~ ^ _1 ). We 
use 2 x 10 5 SPH particles to represent the disc. This provides 
more than enough resolution to capture fragmentation, since 
the Jeans condition (e.g. Bate & Burkert 1997) is easily satis- 
fied everywhere (by more than a factor of one hundred in mass). 
The central star is represented by a sink with radius 2 AU. An 
SPH particle is accreted by the sink if it is within the sink radius 
and is bound to the sink. The central star is allowed to move. 

In both cases the disc relaxes from the initial conditions 
to a quasi- steady state. In Fig. [T] we plot the evolution of the 
disc internal energy with time for both cases. With a uniform 
background radiation field at T ~ 10 K, the disc relaxes to a 
lower- temperature quasi- steady state (bottom line) than in the 
case where radiation from the central star is taken into account 
(top line). Otherwise the evolution of the thermal energy ap- 
pears similar in the two cases. In the following subsections we 
describe each simulation in detail. 

4.1. Disc evolution in a 10 K blackbody background 
radiation field 

In this case, the disc relaxes to a quasi-steady state within 
~ 250 yr; this is the asymptotic phase defined by Boley et al. 
(2006). Thereafter, the disc cools slowly. In Figs. [2]- 14] we plot 
the surface density of the disc, the midplane temperature, the 
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Fig. 1. Internal energy (i) for a disc irradiated by a uniform 10 K 
background radiation field (bottom, red line) and (ii) for a disc 
irradiated by its central star (T BGR (R) ~ R~ l ; top, blue line). 



Toomre parameter g, and the net cooling time, at t = 2000 yr. 
Weak spiral arms form in the disc but they show no tendency 
to fragment (see Fig. [3, despite the fact that the disc is Toomre 
unstable at ~ 30 AU (Fig. [4]). This is because the disc cannot 
cool fast enough; the net cooling time throughout the disc is 
generally > 2t ORB (Fig. [5]), where t ORB is the local orbital period@ 

In Figs. [6] and we plot the azimuthally averaged surface 
density, temperature, Toomre parameter and cooling time (in 
units of the local orbital period) at five times during the disc 
evolution (t = 500, 1000, . . ., 2500 yr). The temperature and 
cooling time are also averaged vertically relative to the disc 
midplane, and the Toomre parameter is calculated using the 
midplane isothermal sound speed. The profiles are essentially 
independent of time, i.e. the disc is in a quasi-steady state. 

This simulation was repeated using a smaller number of 
SPH particles (5 x 10 4 ) and the evolution of the thermal en- 
ergy, surface density, temperature, Toomre parameter and cool- 
ing time were essentially unchanged, indicating that the simu- 
lation is converged. 

4.1 .1 . Comparison with the simulations of Boley et al. 
(2006) and Cai et al. (2007) 

The results of the simulations presented here are similar to 
those of Boley et al. (2006) and Cai et al. (2007). In this sub- 
section, we make a detailed comparison. 

(i) The temperatures we obtain are generally larger than 
those obtained by Boley et al. (2006), on average by ~ 15%. At 
small radii our higher temperatures are probably due to higher 
viscous heating. At large radii (> 30 AU) the temperature differ- 
ences are due to using different background radiation fields: in 



2 In the inner parts of the disc, the ratio of specific heats decreases 
from y 5/3 to y 7/5, due to the increasing temperature and 
the resulting excitation of the rotational degrees of freedom of H 2 . 
Consequently the maximum value of t COOL for fragmentation increases 
to ~ 2t 0RB (e.g. Johnson & Gammie 2003; Rice et al. 2005). 
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Fig. 2. Logarithmic column density in g cm -2 , projected on the 
xy-plane (top), and volume density in g cm -3 on the xz-plane 
(bottom) for a disc irradiated by a 10 K background radiation 
field. 
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Fig. 3. Logarithmic temperature on the xy-plane (top), and on 
the xz-plane (bottom) for a disc irradiated by a 10 K back- 
ground radiation field. 



BGR - 10 K, whereas in the Boley et al. simu- 



our simulation T B( 
lation the disc can cool down to 3 K. 

(ii) In the quasi-steady state, the Toomre parameter in our 
simulation is generally higher than that of Boley et al. (2006). 
This is because in the Boley et al. simulation the temperatures 
are lower, and - due to the strong burst phase, i.e. a phase 
where violent gravitational instabilities occur and rapidly redis- 
tribute angular momentum within the disc - the surface densi- 
ties are higher. However, the region around 30 AU does not fol- 




Fig. 4. Toomre parameter for a disc irradiated by a 10 K back- 
ground radiation field. 




Fig. 5. Net cooling time (in units of the local orbital period, and 
averaged vertically) for a disc irradiated by a 10 K background 
radiation field. 

low this trend, with the Toomre parameter in our simulations 
being smaller (-0.8 compared with ~ 1.45 in Boley et al.); this 
is due to our disc having higher surface density at this radius 
(iii) Our disc relaxes rapidly to its equilibrium state, within 
one disc rotation, without the 'burst' phase reported by Boley 
et al. (2006). We attribute this to two effects. Firstly, the 
higher temperatures in our disc damp gravitiational instabilities 
more rapidly. Secondly, the larger initial density fluctuations 
(oc N~ l/2 ~ 0.0014, where N is the number of SPH particles) 
allow angular momentum to be redistributed by gravitational 
torques without first having to wait for substructure to develop 
by gravitational instabilities. In the Boley et al. simulations, the 
Toomre parameter becomes much lower, and so more violent 
gravitational instabilities develop, and the resulting spiral arms 
re-distribute angular momentum. Cai et al. (2007) and Boley 



3 Note also that Boley et al. use the adiabatic sound speed to calcu- 
late Q, and hence their g-values are increased by a factor y 1/2 ~ 1.3 
relative to ours. 
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Fig. 6. Surface density (azimuthally averaged) and tempera- 
ture (azimuthally and vertically averaged) at five times (t = 
500, 1000, . . .,2500 yr; green, red, blue, cyan, black) during 
the evolution of a disc irradiated by a 10 K background radia- 
tion field. The dotted lines correspond to the initial conditions. 

(private communication) confirm that more powerful ambient 
heating and/or greater noise tend to weaken the burst phase. 

(iv) In the quasi-steady state, our cooling times are sim- 
ilar to those of Boley et al. in the region around 30 AU, but 
somewhat lower outside this region. These differences are at- 
tributable to the differences in surface density and temperature 
noted above. 

(v) The spiral arms that develop in our disc are weaker 
than those reported by Boley et al. (2006). To quantify this 
difference, we decompose the disc structure into a sum of 
Fourier components. We use as our basis a logarithmic spiral, 
R = Roe~ m ^^, where m is the mode of the perturbation, is the 
azimuthal angle of the SPH particle, and f = -m/ tan(J3) is a 
parameter that represents the pitch angle J3 of the spiral (Sleath 
& Alexander 1996). The Fourier transform is then 

I Z! i 6 { u - ln i R j]) $(</>-</>])} e- i( ^ +m(f)) dud(P 

oo %J — ji ■ i 



7=1 



1 N 

- y e- 4 

N 4-f 

7=1 



(an[Rj]+m(f>j) 



(10) 



where (Rj, (pj) are the co-ordinates of particle j. In Fig. [8] we 
plot F(f , m) against f for the m = 1, 2, 3 and 4 modes, for our 
simulation and for the Boley et al. (2006) simulation. For the 
Boley et al. data the sum in Eq. (TTOb is over all points of their 
cylindrical grid, weighted according to the mass mj of each 
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Fig. 7. As, Fig.0 the Toomre parameter (azimuthally averaged) 
and net cooling time (in units of the local orbital period, also 
azimuthally averaged). 

grid cell, i.e. F(£m) = Z%i ™j e<^ R ^ +m ^) I J$ =1 mj. The 
maxima in F(f , m) identify the dominant pitch angles of the 
spirals. In both cases, F has been averaged over ~ 200 yr 
during the quasi-steady state, and in the Boley et al. case F has 
also been divided by 10 to make comparison easier. The spiral 
arms in the Boley et al. disc are about six times stronger than 
in our disc. 

(vi) In our simulation, the central star is represented by an 
accreting sink particle with radius R slUK = 2 AU, which keeps 
the density low near this radius. Consequently we do not see the 
dense rings that Boley et al. (2006) report in the inner regions 
of their disc. 

(vii) In both simulations there is no tendency to fragment. 
Our simulation is also comparable with the simulations by 

Cai et al. (2007) using background radiation fields of 15 K and 
25 K. We estimate that across the disc our temperatures are 
less than 3% higher than in their 15 K case, and approximately 
the same as in their 25 K case. Their discs also show no ten- 
dency to fragment. We conclude that the Indiana University 
Hydrodynamics Group code (e.g. Pickett et al. 1998; Mejia et 
al. 2005) and our SPH-RT code produce very similar results. 
This is significant, because the treatments of radiative transfer 
and hydrodynamics are completely different in the two codes. 

4.1 .2. The vertical temperature profile of the disc 

Our radiative hydrodynamic computational scheme has already 
been tested for the case of collapsing clouds and shown to per- 
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Fig. 8. The mean strengths, F, of the m = 1, 2, 3 and 4 spiral 
modes during the quasi-steady phase, against f = -mj tan(/3), 
where j3 is the pitch angle, for a disc irradiated by a 10 K black- 
body background radiation field (solid lines), a disc irradiated 
by the central star (xlO; dashed lines), and the Boley et al. 
(2006) disc which is irradiated by a 3 K blackbody background 
radiation field (xO.l; dotted lines). 

form well (Stamatellos et al. 2007). Here, we compare our re- 
sults with the analytic prescription of Hubeny (1990) (cf. Boley 
et al. 2007b). 

Hubeny (1990) has calculated the thermal structure of a 
cylindrically symmetric, stationary Keplerian disc. Radiation 
transport is allowed only in the vertical direction, i.e. perpen- 
dicular to the disc mid-plane. Energy is deposited by the shear 
of the Keplerian motion, with a viscosity which is indepen- 
dent of z. Self-gravity, convection and external irradiation are 
all neglected. At any radius, the temperature T at optical depth 
r from the surface of the disc then approximates to 



T(t) = T^<- 



1 

(t - r e ) + — + 



V3 3* r (t)Z 



1/4 



(ID 



where o~ SB T* s is the vertical flux at the disc surface, tq = 
2 k r (o-') cr'dcr' /Z, and cr(z) = p dz. Eqn. (fTTb is derived 
from Hubeny 's Eqn. (8) by setting y j = y H = 1 and w(m) = w. 

We compare the Hubeny solution with our disc at t = 2000 
yr, during the quasi- steady phase. Our simulation accounts for 
external irradiation and for self-gravity, which are not included 



in the Hubeny solution; there are also weak spiral arms in our 
disk. In Fig.|9]the dots represent the vertical temperature struc- 
ture in our disc at different radii (6, 9, 1 1 AU, top plot; 14, 20, 
25, 30, 33 AU, bottom plot); here, temperature is plotted as a 
function of the vertical optical depth from the surface of the 
disc, T(t). The solid lines represent the Hubeny solution at the 
same radii, when T q q is estimated by summing the radiative 
losses L from a thin annulus with cross-sectional area A at each 
radius, and then putting L = 2 A <tsb the resulting Hubeny 
temperatures are lower than in our simulation, especially near 
the central star. However, the Hubeny solution assumes plane- 
parallel symmetry, with a purely vertical temperature gradient. 
In our disc there is also a significant radial temperature gradi- 
ent (see Fig. [6]), and therefore energy diffuses both vertically 
and radially.This results in temperatures near the midplane that 
are higher than the Hubeny solution predicts. To correct the 
Hubeny solution for this effect, we calculate the vertical opti- 
cal depth, r, to the disc surface at each position of the disc, and 
assume that this optical depth also represents how far the radi- 
ation has diffused radially in the disc. We then assume that the 
temperature given by the Hubeny solution at (r, tr - t), where 
tr is the radial optical depth from the inner edge of the disc, ac- 
tually represents the disc temperature at (r, tr). The corrected 
Hubeny temperature profiles (dashed lines) are better fits to our 
disc, apart from the inner inner parts (< 8 AU), where our tem- 
peratures remain higher because they reflect additional heating 
sources that are not included in the Hubeny solution. 

4.2. Disc evolution taking account of radiation from the 
central star 

In reality the disc is likely to be heated by radiation from the 
central star, and this will influence its evolution. In order to 
treat this radiation properly, it is necessary to solve a compli- 
cated transfer problem, which involves both radiation which 
arrives directly from the central star, and radiation which first 
interacts with the more diffuse envelope above and below the 
disc, and is then scattered, or absorbed and re-radiated, towards 
the disc. If a static dust distribution is specified, and radiative 
equilibrium is assumed, this problem can be solved using mod- 
ern Monte Carlo methods (e.g. Wood et al. 2002; Whitney et 
al. 2003; Stamatellos et al. 2005; Pinte et al. 2006; see review 
by Watson et al. 2007). However, if an evolving dust distribu- 
tion is involved, a consistent solution for the dynamics and the 
radiation transport is only feasible - with current computing re- 
sources - if simplifying assumptions are made (e.g. Dullemond 
et al. 2007, and references therein). 

Here we simply assume that radiation from the central star 
can be represented by a blackbody background radiation field 
whose temperature, T BGR , decreases with distance, R, from the 
star. Observations indicate that the disc temperature drops ra- 
dially as R~ q with 0.4 < q < 0.8 (e.g. Beckwith et al. 1990; 
Osterloh & Beckwith 1995). Hence, we set 



T bgr( R ) - 



'■(w) 



+ T 2 



1/2 



(12) 



in Eqn. ©, with T = 1200 K and = 10 K. We refer to this 
prescription as 'implicit stellar irradiation'. (For simplicity, we 
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adopt the same profile for the initial temperature of the matter 
in the disc, cf. Eqn.0. 

The results of this simulation are shown in Figs. I10H15I In 
this case the the disc is initially close to equilibrium, and so its 
temperature hardly changes (see Fig. [14]). Due to the implicit 
stellar irradiation, it is hotter - than the disc irradiated by a 10 K 
background radiation field - and cools fast enough to fragment 
at large radii > 30 AU (Fig. [15]). However, it does not fragment 
because it is not gravitationally unstable (Q > 1.5; Fig.[T5l). 

Implicit stellar irradiation stabilizes the disc. Fourier anal- 
ysis reveals the presence of very weak spiral arms (see Fig. [8]), 
ten times weaker than for the disc irradiated by a 10 K back- 
ground radiation field. This is consistent with simulations of 
discs heated by "envelope-radiation" (e.g. Boss 2001, 2002; 
Cai et al. 2007), and with analytical predictions (Matzner & 
Levin 2005; Rafikov 2005; Whitworth & Stamatellos 2006). 

Comparing our simulation with the T BGR = 50 K case of 
Cai et al. (2007), we find similar g-values across their disc 
and ours. However, our disc is less extended than theirs, since 
they start this simulation from a late phase of their T BGR = 25 K 
case, i.e. after the disc has spread out. We estimate that the m = 
1 and 4 modes are about five times stronger in their simulation 
than in ours (see Fig. [8]). 

5. Summary and conclusions 

We have performed radiative hydrodynamic simulations of the 
inner regions of circumstellar discs with parameters broadly 
similar to those used by Boley et al. (2006) and Cai et al. 
(2007), i.e. a 0.07 M disc around a 0.5 M Q star. Initially, the 
disc extends from 2 to 40 AU, with surface density 2 ~ R~ l/2 , 
and temperature T ~ R~ l . 

We use a new method to treat the energy equation, which 
includes excitation of the rotational degrees of freedom of 
molecular hydrogen, and radiative transfer using realistic dust 
opacities. Our simulations do not have sufficient resolution 
to capture convective energy transport. However, since proto- 
fragments must condense out on a dynamical scale, they do not 
have sufficient time to cool by convection, because this would 
require convective cells to migrate and disperse supersonically 
(Whitworth et al. 2007). Moreover, the surface densities of our 
discs never reach the high values which - according to Mayer 
& Gawryszczak (2007) - lead to proto-fragments cooled by 
convection. 

We have simulated the evolution of a disc irradiated by a 
cool background radiation field (T BGR = 10 K). Despite the fact 
that we use completely different treatments for the hydrody- 
namics and the radiative transport, our results are similar to 
those of Boley et al. (2006) and Cai et al. (2007). Our disc is 
gravitationally unstable at ~ 30 AU, and develops weak spiral 
arms, but it shows no tendency to fragment, because it cannot 
cool fast enough. 

We have also simulated the evolution of a disc taking ac- 
count of radiation from the central star (r BGR ~ R~ l ). In this 
case the disc can cool fast enough to fragment, because of 
its higher temperature, but it does not fragment, because it is 
not gravitationally unstable; it does not even develop notice- 
able spiral arms. This result agrees with previous simulations 



of discs heated by envelope radiation (e.g. Boss 2001, 2002; 
Cai et al. 2007). It also corroborates the analytic predictions 
of Matzner & Levin (2005), Rafikov (2005) and Whitworth & 
Stamatellos (2006). 

Whitworth & Stamatellos (2006) have argued that a mas- 
sive disc will fragment at larger radii (R Z 100 AU), producing 
brown dwarfs and occasionally low-mass hydrogen-burning 
stars or planetary-mass objects. These predictions are corrob- 
orated by the numerical simulations of Stamatellos, Hubber & 
Whitworth (2007). However, if planetary mass objects form at 
such large radii, they are unlikely to migrate inwards to become 
hot Jupiters. More likely they are ejected into the field though 
3 -body interactions. 

We conclude that observed gas giant planets in close orbits 
are unlikely to have formed by gravitational instability. 
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Fig. 9. Vertical temperature profiles, T(t), at different radii (6, 
9, 11 AU, top to bottom, top plot; 14, 20, 25, 30, 33 AU, top 
to bottom, bottom plot). Points correspond to the results of our 
model (azimuthally averaged), the solid lines to the Hubeny 
(1990) solution, and the dashed lines to the Hubeny solution 
corrected for the fact that radiation diffuses both vertically and 
radially. 
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Fig. 12. As Fig. HI but taking account of radiation from the cen- 
tral star. 



Fig. 10. As Fig 
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Fig. 13. As Fig. [3 but taking account of radiation from the cen- 
tral star. 



Fig. 11. As Fig.0 but taking account of radiation from the cen- 
tral star. 
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Fig. 14. As Fig.0 but taking account of radiation from the cen- 
tral star. 




Fig. 15. As Fig. [71 but taking account of radiation from the cen- 
tral star. 



